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ABSTRACT 

In  this  paper,  a  new  software  availability/reliability  model  is  developed  where  life¬ 
times  and  repair  times  have  general  system-state-dependent  distributions.  Multiple  errors 
may  be  introduced  or  removed  through  repairs.  The  model  is  formulated  as  a  multivari¬ 
ate  Markov  process  and  contains  many  other  models  appeared  in  the  literature  as  special 
cases.  The  exponentiality  assumption  prevalent  in  the  literature  is  totally  eliminated. 
Expressions  of  various  performane  measures  of  practical  interest  combining  availability 
and  reliability  of  the  software  system  at  time  t  are  derived.  Using  the  matrix  Laguerre 
transform  of  Sumita(1984),  corre  ponding  computational  procedures  are  also  developed. 
A  numerical  example  is  given,  demonstrating  speed,  accuracy  and  stability  of  these  pro¬ 
cedures. 
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§.0  Introduction 


The  price  performance  revolution  of  computer  hardware  has  been  dramatic .  w  t 
cost  of  labor  has  been  steadily  increasing.  Consequently  the  production  and  maintcnam  ■ 
cost  of  software,  in  contrast  to  that  of  hardware,  has  been  rapidly  growing  and  has  become 
one  of  central  issues  in  system  design.  Dating  back  to  the  late  60’s  substantial  research 
efforts  have  been  devoted  to  the  study  of  software  failure  phenomenon  and  the  prediction 
of  software  performance.  Two  recent  survey  papers  by  Ramamoothy  and  Bastani(1982) 
and  Shanthikumar(l983)  contain  approximately  150  refereneces  on  the  issues. 

Until  recently,  however,  the  effect  of  multiple  error  generation  and  removal  from  the 
system  during  the  repair  has  not  been  properly  incorporated  in  the  literature.  Kre- 
mer(1983)  has  derived  a  software  reliability  model  where  the  number  of  software  system 
increases  or  decreases  by  at  most  one  during  repairs.  He  has  provided  performance  mea¬ 
sures  for  this  model  using  the  results  available  for  non-homogeneous  birth-death  process. 
Sumita  and  Shanthikumar(l984)  have  developed  a  general  Markov  chain  model  where 
multiple  errors  may  be  introduced  or  removed  from  the  system  during  repairs.  Assuming 
that  the  software  failure  rate  is  propotiona!  to  the  number  of  software  errors  present  in 
the  system,  expressions  for  various  software  reliability  measures  of  interest  are  derived 
and  corresponding  computational  procedures  are  developed. 

The  exponentiality  assumption  employed  in  the  model  of  Sumita  and  Shanthiku- 
mar(  198-1)  is  rather  restrictive.  To  preserve  the  Markov  chain  property,  for  example, 
the  software  repair  time  is  assumed  to  be  negligible,  ignoring  the  availability  of  software 
system  at  time  t.  The  problem  of  a  combined  availability /reliability  analysis  of  software 
system  was  addressed  by  Shant hikumarj  1984)  in  a  limited  model  where  the  number  of 
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software  errors  is  either  unchanged  or  reduced  by  one  during  the  repairs.  The  purpose 
of  this  paper  is  to  develop  a  general  multivariate  Markov  mode!  for  a  software  system 
with  multiple  software  error  generation  and  removal  during  the  repairs.  The  exponential¬ 
ly  assumption  is  totally  eliminated  and  general  repair  time  distributions  are  explicitly 
incorporated.  E’  .essions  of  various  performance  measures  combining  availability  and 
reliability  of  s<-  cware  system  at  time  t  are  derived.  Using  the  matrix  Laguerre  transform 
of  Sumita(x984)  as  a  key  tool,  corresponding  computational  procedures  are  also  devel¬ 
oped.  Many  other  models  appeared  previously  in  the  literature  can  be  treated  as  special 
cases  of  this  model. 

In  Section  1,  we  develop  a  new  software  availability /reliability  model  having  system- 
state-dependent  lifetimes  and  repair  times.  Multiple  errors  may  be  introduced  or  removed 
through  repairs.  The  model  is  formulated  as  a  multivariate  Markov  process  and  does  not 
require  exponentiality  at  all.  By  studying  the  probabilistic  flow  in  the  corresponding  state 
space,  various  time  dependent  entities  are  analyzed.  In  Section  2,  expessions  of  man} 
performance  measures  are  derived  in  terms  of  these  probabilistic  entities.  Performance 
measures  combine  time  dependent  availability  and  reliability  of  the  software  system. 
Computational  procedures  for  evaluating  these  performance  measures  are  developed  in 
Section  3,  using  the  matrix  Laguerre  transform  of  Sumita(1984).  Section  4  is  devoted 
to  numerical  implementation  of  the  procedures  demonstrating  their  speed,  accuracy  and 

stability. 
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We  consider  a  software  system  which  contains  several  software  errors.  These  software 


errors  cause  software  system  failure  time  to  time.  Upon  failure,  the  software  system  is 
repaired.  During  repairs,  multiple  software  errors  may  be  introduced  or  removed.  More 
formally,  the  following  assumptions  are  incorporated  in  our  model. 

(ASl)  The  maximum  number  of  software  errors  in  the  software  system  is  limited  to 
K,  0  <  K  <  oc.(If  the  maximum  number  of  errors  exceeds  A,  then  the  performance  of 
software  system  becomes  untolerable  and  the  system  would  be  discarded.) 

(AS2)  At  time  t  =  0,  the  software  system  starts  functioning  and  there  are  A'  errors  in 
the  software.  Here  Ar  is  a  discrete  nonnegative  random  variable  with  probability  vector 
bT  =  ( 60 , . . .  ,6k)  where  6,  =  P\ X  =  i],  i  =  0, 1, . . . ,  K. 

(AS3)  If  there  are  n  software  errors  upon  completion  of  a  repair,  then  the  time  until 
next  software  system  failure  has  c.d.f.  An(x),  n  =  0, 1, 2, . . . ,  K .  In  particular  -4q(x)  =  0, 
x  >  0,  i.e.  if  there  is  no  error  in  the  software  system,  then  there  will  be  no  software 
system  failures.  It  is  assumed  that,  for  n  >  1,  An(x)  is  absolutely  continuous  with  p.d.f. 
an(x)  and  hazard  function  r?n(x)  =  an(x),' A„(x)  where  An(x)  =  1  -  A„(x). 

(AS4)  If  there  are  n  errors  at  the  begining  of  a  repair,  then  the  repair  time  has  a 

c.d.f.  Rn(x),  n  =  1,2 . A'.  It  is  assumed  that  An(x)  is  absolutely  continuous  with 

p.d.f.  rn(x)  and  hazard  function  ?„(x)  =  rn(x)/ A„(x). 

(AS5)  The  probability  that  there  are  n  errors  remaining  in  the  software  system  im¬ 
mediately  after  a  repair  given  that  there  were  k  errors  in  the  software  system  just  before 
the  begining  of  the  repair  is  p^n.  For  notational  convenience,  we  define  poo  -  1,  Pon  =  0, 

n  >  I. 
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(AS6)  All  software  system  lifetimes  and  repair  times  are  mutually  independent. 

This  model  can  be  expressed  as  a  multivariate  Markov  process  in  terms  of  the  follow¬ 
ing  stochastic  processes: 

!0,  if  the  software  system  is  under  repair  at  time  t. 

1,  if  the  software  system  is  functioning  at  time  t. 

(1.2)  A/(f)  =  the  number  of  software  system  failures  occured  in  [o,t)- 

(1.3)  iV(t)  =  the  number  of  errors  in  software  system  at  time  t. 

(1.4)  A'(t)  =  the  elapsed  time  since  the  last  transition  of  /(f),  i.e., 

X(f)  =  t  -  r,  where  r  =  sup  {x  :  |/(x  +  )  -  /(x-)|  =  1} 

0<x  <t 

From  (AS2),  one  has  A'(0~)  =  0.  For  notational  convenience  we  assume  that  J(O-)  =  0 
and  7(0+)  =  1.  Clearly  the  multivariate  process  '  f(f),  A/(f) ,  N(t) ,  A”(f )  j  is  Markov.  The 
state  space  of  the  multivariate  process  and  its  typical  transition  behavior  are  depicted  in 
Figure  1.1. 


/(0  =  0  I(t)  =  1 


Figure  1.1 


Let 


F.m.n(x,0  =  P[  I(t)  =  z ,  J \f(t )  =  m.  Ar(0  =  n.  -Y(0  <  xj, 


and  define 


(1.6)  ~  (•r'  0  ■ 

The  partial  differentiability  /imn(x,t)  =  ^ Ft,m,n  (x,  t)  can  be  shown  through  a  renewal 
argument  (see,  e.g.,  Qinlar(l969)),  except  the  case  m  =  0  for  which  generalized  densities 
/l ,0,n (x,  t)  =  bn6(t  —  x)An(x)  are  involved.  As  we  will  see  in  Section  2,  all  performance 
measures  of  practical  interest  for  the  study  of  availability/  reliability  of  the  software 
system  can  be  expressed  in  terms  of  (1.5)  and  (1.6).  In  the  remainder  of  this  section, 
we  derive  transform  results  of  (1.6)  by  applying  the  state  space  method  of  Keilson  and 
Kooharian(1960,  1962).  The  corresponding  computational  procedures  will  be  discussed 
in  Section  3. 

We  observe  from  Figure  1.1  that  for  the  process  to  be  at  (z,m,n,x)  at  time  t  where 
0  <  i  <  t,  the  process  must  have  entered  (i,m,n,0)  at  time  t  —  x  and  has  remained 
in  the  state  (z,  m.n,  )  for  the  length  of  x.  We  note  from  (AS3)  that  /lo(x)  —  0,  x  >  0. 
Therefore  once  the  process  enters  (l,m,0,  ■),  it  remains  there.  Hence  one  has 


fo,m,n  (x,  f )  -  /o.m,n(0~,f  -  x)/?n(x),  m,  Tl  >  1 


-  /l,m,n(0  + ,  t  -  l)/lR(x),  m,Tl  >  0. 


Bv  a  similar  argument,  boundary  conditions  can  be  found  as 


/, ,m ,rt(x,0)  =  0,  i  =  0,  1,  m,  n  >  0,  x  >  0, 
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(1.10) 


/o.m.n(0-,0  =  J  /i,«-i,«(x,0»7n(x)dx.  m,n  >  1. 


Futhermore  one  has 


(1.11) 


/l.m.«(0  +  ,0  =  E*  Pknfofo.m.k(x'tUk(x)dx'  m  >  l,n  >  0, 

/l,O.n(°+»0  =M(0^ 


Substituting  (1.7)  and  (1.8)  into  (1.10)  and  (1.11),  one  obtains  for  t  >  0, 

(1.12)  .  /o,m,n(0-,0  =  J‘  -  x)an(x)dx,  m>l,  M>1, 


(1.13) 


/i,m,n(0-,0  =  E*  P*n/o/o.m.*(0-,t  -  x)r*(x)dx,  m  >  l,n  >  o 

/l,0,n(0^,t)  =  M(0,  n>°- 


(  £»,m,n(0  — ,  s)  —  /o’  C  "Vt,m,n(^Td)^ 

Equation  (1.12)  and  (1.13)  can  be  best  described  in  terms  of  the  transfrom  of  (1.14)  using 


the  matrix  notation.  For  notational  convenience,  we  introduce 


(1.15) 


/o,m,o(x,<)  -  0,  m  >  0,  x,t  >  0 


and  define  the  transform  vectors 


(1.16) 


^m(0 i,i)  =  [^,.m1o(0-i',;)>---  i£i,m,K'(0-,s)  ] 

6  (iy,s)  =  (^,,m.o(u’’5)’-"  ’  iKt.m.K^5)  ]• 


Let  diag{co, . . . ,  ck)  be  a  (K-f  1)  x  (K  +  1)  diagonal  matrix  whose  n-th  diagonal  element  is 
cn- 1-  We  define  QD{s)  =  diag{  0,  qi  (s) . (s)  }  and  £D(s)  --=  dia9{0,  Pl  (s), . . .  ,  PK  (*) 


4  JV.  .  JV  -.-V  .  ■V  ■- 


-vV  ,  „*\  > 


v  -s*—  tj-  ■"T'Tr'  "■ 


where  aril's)  =  /0°°  e  3zan(x)dx  and  pn(s)  =  /0°°  e  3Irn(x)dx.  One  then  sees  from  (1.12) 
and  (1.13)  that 


(1.17) 
and 

(1.18) 


£jm(°  ~5)  =  ^ 


±lj0  +  'S)  =^.m(°-’S)£D(S)£’ 

^f,0(0+’5)  =  ^ ' 


m  >  1, 


where  P  =  (p,,).  Hence  one  concludes  that 


(1.19) 

Vi 

.S'.'r 

k*  .> ; 

and 

►7  A 

K«>: 

(1.20) 

(£?  (0-r 

— l,m ' 

^  {^D^S)PD^^)m  *  =d('S)’ 


(1-20)  £[,m(°  +  ’5)  =  ^ 

Finally  from  (1.7),  (1.8),  (1.19)  and  (1.20).  one  has 

(L21)  =  J^M.rn(0~'S)U-  PD{*  +  «)), 

and 


(1.22) 


-  T  1  r 


=  7~ - ;  £i,m(°~’S^  =  _  £/?(*"  *  5)),  m  >  0. 


We  note  that  the  spectral  radius  of  y£_D[s)p  [s)Pj  is  strictly  less  than  one  for 
Re{s)  >  0  and  £*=0  (qd($V  (s)£)m  ~  (l  ~  0D(s)p  (s) P)  '  .  Hence 


(1.23) 


E  £l.m(0*’5)  =  hT \l  -  ^o('s)^(5)£i  ‘od(sK 


(1.24) 


E  £[,m(0^5)  =  6r:/-aD(S)£  (S)£)- 


1-ViAi  r_',  J»_*  xA  .V.*,  iv.1  .V.’.  M 


Let  1  be  a  vector  of  length  (A'  +  1)  having  all  elements  equal  to  one.  It  can  be  readily 
seen  from  (1.21)  through  (1-2-1)  that 


This  model  contains  the  model  of  Sumita  and  Shanthikumar  (1984)  as  a  special  case 
where  .4n(x)  =  1  —  e~nXz  and  An(x)  =  1  for  x  >  0.  Hence  many  other  models  contained  in 
Sumita  and  Shanthikumar(1984)  are  also  contained  in  this  model.  Shanthikumar(l984) 
has  discussed  another  special  case  of  this  model  where  pt;  >  0  if  and  only  if  j  —  i  or 
3  =  *“  1- 
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§.2  Performance  of  the  Software  System:  Availability  Reliability  Measures 


Ln  this  section  vve  introduce  various  availability  reliability  measures  for  the  study  o! 
performance  of  the  software  system.  Based  on  the  analytical  results  of  Section  1,  we 
derive  explicitly  expressions  of  these  availability/reliability  measures. 

A  natural  availability;  reliability  measure  is  the  joint  probability  of  I(t)  (whether  the 
software  system  is  working  or  not  at  time  t),  \f(t)  (the  number  of  software  system  failures 
occured  in  !0,  t)  ),  and  X{t)  (the  number  of  errors  in  the  software  system  at  time  f).  One 


P  I{t)  =  t,  M{t)  =  m,  X(t)  =  n  1  =  Ft,m,n(- ocj). 


The  joint  probability  of  I(t)  and  X(t)  is  also  of  interest.  For  notational  convenience,  we 
define  iv),o,n(~oc,  t)  =0,  n  >  0,  t  >  0.  One  then  sees  that 


P'I(t)  =  *,  X(t)  =n j  =  £  Ftifn,n(  +  oc,0. 


The  availabilitv(unavailabilitv')  measure  is  the  probability  that  the  software  svstem  is 
functioningfnot  functioning)  at  time  t,  given  by 

K  co 

(2.3)  P\I{ t)  =  i]  =  Yi  Fx,m.n(-° C'O- 

n  =  0  m  =  0 

The  time  until  complete  debugging  is  an  important  performance  measure.  Let  .V Vl 
be  the  time  required  for  completely  eliminating  all  errors  in  the  software  system  giver, 
that  Ar(0  — )  =  k  with  probability  6*.  We  denote  the  distribution  function  of  A  y,  i<\ 
5v(o)o(r)  =  P\ Xs(o)o  -  x\-  Since  the  states  (l,m,0,  ■)  are  absorbing,  one  finds  that 

OC 

(2.4)  5,v(o).o(0  =  PdV(0  =  °i  =  Y.  Pi,m,u(-x,t). 
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In  the  study  of  performance  of  software  systems  it  is  of  more  practical  interest  to 
predict  the  software  availability/reliability  based  on  the  past  observation  of  the  system. 
In  order  to  obtain  performance  measures  in  this  context,  we  impose  the  following  history 
of  the  system: 

(CD)  A  repair  has  just  been  completed  at  time  to  and  there 
have  been  m  software  failures  in  jO,  £ o ] - 

The  probability  distribution  of  the  number  of  errors  remaining  in  the  software  system 
at  time  to  given  (CD)  plays  a  crucial  role  in  this  analysis.  Let  ffkitutTn  =  P[N(to)  = 
k\{CD)]  and  define  . . .  ,0K\uum  !•  0ne  then  finds  thal 

(2-5)  = 

The  vector  &Tti,  m  fully  describes  the  state  of  the  system  at  time  t0  under  (CD),  which 
then  provides  an  initial  state  probability  vector  for  the  system  behavior  after  time  to-  To 
emphasize  the  dependence  of  F,im  n(i,t)  on  the  initial  distribution  6r  we  write 

(2.6)  Ft,m<n(x,t\b)  = 


P\I(t)  =  :,  M(t)  —  rn,  .V(f)  =  n,  A'(t)  <  zr|Ar(0-)  =  k  with  probability  6*]. 

Given  (CD),  the  joint  probability  of  the  system  availability  /unavailability  and  the  num¬ 
ber  of  errors  in  the  software  svstem  at  time  tn  -f  r  is  then  obtained  as 


P{I(t0  t  r)  =  :,  Ar(t0  tr)  =  n|(CD)]  =  V  FtiJ,n(~oc,  rj  f  (  J. 

j=o 


Similarly  the  joint  probability  of  the  system  availability/unavailability  and  the  number 
of  software  failures  that  may  occur  in  [(o,*o  -r  r)  under  (CD)  is  given  by 


(2.8)  P[I(t0  -r  r)  =  i,  M(to  -f  t)  -  m  =  j\(CD)\  =  FtiJ<n(^oo,  r|  0  ). 
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Let  Y'v({  )  o  denote  the  random  duration  required  for  completely  eliminating  the  software 
errors  given  (CD).  The  corresponding  distribution  function  is  defined  by  G  sr[to),oiT)  = 
P\YSVl>lo  5:  T'] ■  ^  can  be  readily  seen  that 

OO 

(2-9)  G,v(fo)i0(r)  =  Y  F,,;io(+oc,r!5ji(i  m). 

j=o 

Finally  an  important  reliability  measure  is  the  time  until  the  next  software  system  failure 
given  (CD).  We  denote  this  random  variable  by  T\to  m.  The  survival  function  U-j,u  m(r)  = 

P\Tilo,m  >  r]  is  then  S>ven  b>' 

__  K  _ 

(2-10)  »'!t0.m(')  =  £&[  to,mAk(r). 

k=0 

We  note  that  T]h>m  is  dishonest  and  W|t()  m(r)  —  /?0( to>m  as  r  —  oo. 

It  should  be  noted  that  the  marginal  process  N(t)  is  absorbing  and  all  the  perfor¬ 
mance  maesures  described  above  are  time  dependent.  One  time  independent  performane 
measure  of  interest  is  the  distribution  of  the  number  of  software  system  failures  that 
occur  before  all  software  errors  are  completely  eliminated.  We  denote  the  correspond¬ 
ing  randon  variable  by  D.  One  then  sees  that  dm  =  P\D  —  ml  =  Fi,m,o(-oc,  -foo)  = 
limf—o-t-  5 Pi,m,o(0> •$).  Hence  one  has 

(2.1 1)  cfm=4r(rZ)%  m  >  0. 

where  V  —  ap  {0,1,1,...,!},  and  ej  -=  (1,0.0. ...  ,0). 
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§.3  Computational  Procedures 

We  have  seen  in  Section  2  that  all  performance  measures  can  be  obtained  if  £,^(+00,  t ) 
and  (0-^,t)  are  computed.  In  this  section,  we  develop  numerical  procedures  for  eval¬ 
uating  these  probabilistic  entities.  We  assume  that  both  an(i)  and  r„(x)  belong  to  the 
class  of  rapidly  decreasing  functions  of  Dym  and  \lckean(1972),  so  that  the  corresponding 
Fourier-Laguerre  coefficients  are  also  rapidly  decreasing,  see,  Keilson  and  Nunn(l979). 
The  Laplace  transform  of  Fjm ( -i-  oc ,  t )  denoted  by  m(s)  =  /0°°  e~stFjm[  +  oo,  t)dt  can 
be  found  from  (1.21)  and  (1.22)  by  setting  tv  =  0-f,  that  is 

io.m('S)  =  (l/s)gm{0^*)iL  -  PD{* )).  ™>l, 

(3<1)  .r 

>  $l,m(s)  =  ^  °- 

The  entities  ~,s)  and  £^(0  —  ,s)  needed  here  are  given  in  (1.19)  and  (1.20).  The 

inversion  of  these  transforms  require  multiple  convolutions  of  matrix  functions  in  the 
time  domain.  The  matrix  Laguerre  transform  developed  by  Sumita(l984)  provides  a 


computational  vehicle  for  this  purpose. 

The  Laguerre  transform,  introduced  in  Keilson  and  Nunn(l979)  and  Keilson, Nunn 
and  Sumita(l98l)  and  further  studied  by  Sumita(198l),  provides  an  algorithmic  frame¬ 
work  for  the  computer  evaluation  of  multiple  convolutions  and  other  continuum  opera¬ 


tions.  The  transform  based  on  generalized  Fourier  series  employs  the  Laguerre  functions 
as  a  basis,  and  maps  the  functions  /(x)  in  L 2  into  discrete  sequences  (/“)!“.  Corre¬ 
spondingly,  various  continuum  operations  are  mapped  into  lattice  operations,  thereby 
providing  the  desired  algorithmic  basis.  Recently  the  formalism  has  been  extended  to 
the  matrix  form  for  the  study  of  semi-Markov  processes  by  Sumita(  1984)  which  is  the 
crucial  numerical  tool  employed  here.  For  the  reader  s  convenience  a  concise  summary  of 
the  matrix  Laguerre  transform  is  given  in  the  Appendix.  The  notation  there  is  employed 
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throughout  the  rest  of  the  paper. 

Let  an(x)  and  rn(x)  have  the  Laguerre  sharp  coefficients  (a^.k)kL0  and  (r*j.)£l0.  One 
then  has 

(3.2)  cn(x)  =  J  an(x  -  y)rn [y)dy  ~  c%k  =  ^  at.k-,rt. ,  • 

;=0 

Let  g(x)  be  the  matrix  function  defined  by 

I- CO 

(3.3)  yo  e“',f£(x)dx  =  £d(s)£d(s)P. 

It  can  be  re\dilv  seen  that  g(x)  has  a  sequence  of  Laguerre  sharp  coefficient  matrices 

I  a  ! 

£*  =  r*-„.  where 

(3.4)  **“„  =  0  <  t,j  <  K,  k  =  0,1,.... 

We  note  that  an(x)  =  ro(x)  =  0  so  that  =  r*k  =  0  for  £  >  0  and  hence  = 
0,  A:  >  0.  If  £!m|(x),  having  the  Laguerre  sharp  coefficint  matrices  (gf  (m)),  corresponds 
to  (£D('s)P£)(s)P)rr\  one  has 

(3.5)  iT(m  +  1)  =  T,  £;*■ 

J=0 

Hence  the  Laguerre  sharp  coefficient  vectors  of  /^"  (0-r,f)  can  be  easily  computed  us¬ 
ing  (1.20).  Those  corresponding  to  0~J)  can  then  be  obtained  by  convolving 

the  resulting  sharp  coefficient  vectors  for  with  the  Laguerre  sharp  coefficient 

matrices  (a^J2°  for  gD(x)  corresponding  to  qD(s).  It  should  be  noted  that  a^k  is  a  diag¬ 
onal  matrices  whose  n-th  diagonal  element  is  a*_,  Once  the  Laguerre  sharp  coefficient 

j* 

vectors  are  found,  the  function  values  of  /lm(0 J)  can  be  computed  straightforwardly 
following  the  inversion  procedure  described  in  the  Appendix. 
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Evaluation  of  the  function  values  of  Fi  m(  +  oc,£)  requires  more  caution  since 
+  oo,t)  may  not  be  in  Li.  In  particular  /']im,o(4-oo,  t)  converges  to  a  positive  constant 
(see(2.11))  as  t  — *  oo.  From  (3.1),  one  sees  that  Fjm(a-oc,t)  is  differentiable  with  respect 
to  t  for  all  m  >  0.  Let 

(3-6)  =  ^£,rm(^oc,0 

and  define 


It  should  be  noted  that  m(s)  =  s  ^irn(s)  -  m(+oo, 0).  Clearly  F^m(  +  oo,0)  =  0r 

for  t  =  1,  m  >  1  or  t  =  0,  m  >  0  and  F1)0(-<-oo,  0)  =  b.  Hence  from  (3.1)  one  obtains 

(3  g  |  ~  £D(S))>  m  ^ 

I  2 f,m(5)  =  £l!m(0  +  ’5)U  “  gtf(S))  “  60mf,  TTl  >  0. 

Here  <5om  =  1  if  m  =  0,  bom  =  0  otherwise.  It  can  be  readily  seen  from  (3.8)  that 

are  Laguerre  transformable  (see  Remark  3.1)  and  the  corresponding  Laguerre 

sharp  coefficient  vector  denoted  by  (g*m.k)% L0  can  be  found  as  before.  Let 

(3-9)  £.rm(0  =  I" 


We  note  from  (1-20)  that 


(3.10) 


=  Fd'gr-'V  <  Ir.  mil, 
?L  =  =  fujgr  <iT.  m>  o. 


Here  the  n-th  conponent  of  m  is  the  probability  that  upon  the  occurence  of  the  m-th 
failure  there  are  n  software  errors  in  the  software  system.  Since  there  may  be  positive 
probability  that  all  software  errors  are  eliminated  before  reaching  the  m-th  failure,  it  is 
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possible  lo  have  qj ^  1  <  1.  The  vector  qj  m  has  a  similar  probabilistic  meaning.  One  now 
sees  from  (3.8)  and  (3.10)  that  since  <?  m.o  —  0, 

|  2T  (0)  =  0,0 . O’,  m  >  1 

(3 - i 1 ) 

(  =  .yi.m,0,O - ,0)  -  SombT.  TTX  >  0. 

___7'  ^  j- 

Hence  Gim(0)  =  2,  m(0)  <  00  an^  Equation  (3.9)  is  well  defined.  By  applying  one  of 
the  operational  properties  of  the  matrix  Laguerre  transform  (see  ( A 1 3 ) ) ,  the  Laguerre 
sharp  coefficient  vectors  of  G,  m(f)  can  be  generated  from  {g“m )£!,-,  for  Oi.e  then 

finally  has 

(3-12)  £?m(-oc.O  =G^(0)-Girm(0-drnU«5,i^- 

Remark  3,1 

We  note  from  (3.8)  that  gj n{t)  is  a  generalized  vector  function  involving  the  delta 
function  <5(x).  Hence  the  components  of  g[0(t)  are  not  in  L;.  In  a  recent  paper  by  Keilson 
and  Sumita(198-l),  it  has  been  shown  that  the  Laguerre  sharp  transform  exists  for  any 
finite  signed  measure  preserving  all  basic  operational  properties.  The  Laguerre  sharp 
coefficients  for  g?  (f)  therefore  exists. 


§.4  Numerical  results 

In  this  section,  we  demonstrate  the  efficiency  of  the  computational  procedures  devel¬ 
oped  in  Section  3  through  a  numerical  example.  Tables  and  graphs  illustrating  numerical 
results  are  given  at  the  end  of  this  section.  We  consider  a  software  system  with  the 
following  features  corresponding  to  (ASl)  through  (AS5)  of  Section  1. 

(4.1)  K  =  4. 


b(  =  (0.2.  0.2,  0.2,  0.2.  0.2) 


Ao(x)  =  0  and  An(x)  =  f  ~  .  -y*  ne  2ydy ,  1  <  n  <  4. 

Jo  (4  -  n  . 


M*)  =  f 

JO 


(4  -  n)! 


y 4_ne_3ydy,  1  <  n  <  4. 


1.000 

0.000 

0.000 

0.000 

0.000  \ 

0.700 

0.150 

0.100 

0.030 

0.020 

0.300 

0.425 

0.150 

0.100 

0.025 

0.125 

0.225 

0.400 

0.150 

0.100 

0.075 

0.125 

0.200 

0.400 

0.200  / 

We  note  that  the  lifetime  and  the  repair  time  of  the  system  when  there  are  n  software 
errors  in  the  system  are  the  sum  of  (5  -  n)  independent  exponential  random  variables 
with  parameter  2  and  3  respectively. 

The  Laguerre  sharp  coefficients  (a*)2°  of  an  exponential  function  an(x)  =  6e~Sl , 
6  >  0.  can  be  found  analytical)-  (see  Keilson  and  Nunn(198l))  as 


ao  = 


„*  = _ L_  LJ 


,  n  >  1. 


This  formula  enables  one  to  generate  the  Laguerre  sharp  coefficient  matrices  correspond¬ 
ing  to  «c({)  and  p  (s)  via  discrete  convolutions.  In  actual  computation,  the  first  101 


1G 


t 


coefficients  were  used.  Using  equation  (1.17)  and  (1.18),  the  Laguerre  sharp  coefficient 
vectors  of  /^(0-S f),  i  =  0.1,  can  then  be  obtained  via  discrete  vector-matrix  convolu¬ 
tions.  As  discussed  in  Appendix,  the  moment  formula  of  the  matrix  Laguerre  transform 
provides  a  huristic  tool  for  checking  accuracy.  Let 

rcc 

(4.7)  £jk)  =  fo  0  <  ft  <  2. 

Byd  ifferentiating  (1.17)  and  (1.18)  with  respect  to  s  at  s  =  0,  one  then  finds  the  following 
recursion  formulas: 

'  =  h- 

(4.8)  ^JO)  =  /fo.JO)rC 

.^(°)  =  f£1,m-1(o)r- 

’i(i)  =  KJl)^i(°)Wi))£ 

(4.9)  1  Mf0(l)  =  Q 

'  *£  J  =  )  =  (^(2)r  -  2pJm(1)A£KFD(1)  -  ffI.m(0)MHPo(2))  £ 

(4.10)  p[0(2)-Q 

.  ffj,rr.(2)  =  ~  “ff  l,m-l  (^=£lFD  (  ^  ~ 

Here  M.LFD^)  =  dla9{0-  fo°  z*ai  (x)dx,  ...,  f0°°  x*a4(x)dx}  for  0  <  /c  <  2  and  .\/ppD(/c) 
is  defined  similarly  for  repair  time  distributions.  Using  (4.8)  through  (4.10),  njm(k)  were 
calculated  for  i =  0,1,  1  <  m  <41  and  0  <  k  <  2.  These  values  were  then  compared 
with  the  values  obtained  from  the  Laguerre  sharp  coefficient  vectors  of  fT  (O-r.f)  and 

— t,m  v  ’  ' 

the  moment  formulaof  (A15).  The  relative  errors  were  found  to  be  bounded  by  1  x  10~‘° 
in  this  range  of  t,  m,  and  k ,  providing  excellent  accuracy.  In  Table  4.1,  this  comparison 
is  exhibited  for  z  =  1  and  m  =  10.  The  value  rn  =  10  will  be  used  subsequently  for 
evaluating  conditional  performance  measures. 


„  J'm 
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The  Laguerre  sharp  coefficient  vectors  of  gTm  (f)  in  (3.6)  can  be  found  from  those  of 
(0  +  ,  f)  with  one  additional  discrete  vector-matrix  convolution  using  (3.8),  which  in 
turn  lead  to  Laguerre  dagger  coefficient  vectors  of  G,,m(0  of  (3.9)  via  the  operational 
property  of  the  matrix  Laguerre  transform  given  in  (A13).  The  values  of  Fjm[+oc,  t ) 
needed  can  then  be  calculated  from  (3.12).  Since  most  of  the  performance  measures 
involve  the  expressions  Dm=i£o,m(ioo«!)  and  Flm(+oc,  t),  we  generate  the  La- 

_ T 

guerre  coefficient  vectors  of  S,,,\/(0  defined  by 

(4-11)  S.vr(0=  £  Cm(0- 

m  =  1  -t 

Both  the  Laguerre  transform  and  the  tail  integral  operation  are  linear  and  this  can  be 
accomplished  by  merely  adding  the  Laguerre  sharp  coefficient  vectors  of  qT  (0-r,t)  over 

— i,m '  ' 

m,  0  <  m  <  M,  and  then  applying  the  operational  property  of  (A  13)  to  the  resulting 
sum. 

To  check  the  accuracy  of  the  truncated  Fourior-Laguerre  transform  representation  of 

_ 7* 

one  again  uses  the  moment  formula  of  (A15).  Let 

(4-12)  dmW  =  }“  >“  iljtW- 


By  differentiating  (3.8)  with  respect  to  s  at  s  =  0,  one  finds  that: 

d.M  =  £jm-n 

dJ°)  =  -  D  -  WT 


(4.13) 


(4.14) 


(4.13) 


Ei t ' )  =  fii™  ( ■ I )  ( L  -  r )  -  <m  (0)  nRPD  ( 1 ) 
\J.r.W=‘ilJm-r)-!LlJo)VLFD(i) 


'X»P)  =  dS-U  -  n  -  *d.jmSPD{  1)  -  d.jmRPDW 

nU->  =  d.S-U  -  £')  -  2ef.„(i)-'Ji„(i)  -  ttfm(0).'ltfP(2). 
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It  should  be  noted  that 


(4.16)  f°°£u(t )dt=  T  rTm(i);  2  rts^(t)dt 

Jo  '  , 

■n=  1  -i 


M 

E  : 

m  =  1  -t 


The  zero-th  and  the  first  moment  of  51|.\/(0  were  calculated  using  the  Laguerre  sharp 
coefficient  vectors  and  the  momemt  formula  of  (A15).  These  values  were  then  compared 
with  values  of  r,>m(l)  and  £m=i-t  L^m(2)  generated  from  the  recursion  formu¬ 

las  in  (4-13)  through  (4.15).  Both  the  relative  and  absolute  errors  decrease  in  XI  and 
those  errors  were  found  to  be  bounded  by  1  x  10" '  and  3  x  10~B  respectively  for  XI  =41. 
This  comparison  is  summerized  in  Table  4.2.  The  truncation  level  M  has  cross  examined 
using  the  sequence 


A/ 

(4.17)  vM{t)  =  E{££m+i<  +  oc,0  +  £f.mK° c,t))l. 

m=u 

which  converges  to  1  as  XI  — *  oc  for  all  t  >  0.  For  the  value  M  =  41.  (1  —  r.\/(t))  was 
found  to  be  bounded  by  1  x  10~18  for  0  <  t  <  20.  as  exhibited  in  Table  4.3. 

The  time  independent  performance  measure,  dm,  of  (3.11)  describing  the  probability 
of  having  m  software  failures  before  complete  debugging  is  depicted  in  Figure  4.4.  In 
Figure  4.5,  the  compound  availability/reliability  measures  P./(f)  —  1,  XI (t)  =  3,  Ar(t)  = 
n]  =  i'i,3,n(  +  oo» t)  are  plotted  for  0  <  n  <  4  and  0  <  t  <  20.  One  observes  that 
limr_oo  Pi ,3,0 (*f- oc,  t)  —  =  0.14259  and  lim(_oo  ^1,3, n(  —  oc.  t )  =  0.  1  <  n  <  4.  as 

expected.  The  joint  probability  P\I(t)  =  t,  A(f)  =  n  ot  (2.2)  are  exhibited  in  Figure 
4.6  and  4.7  for  t  =  0,  1  <  n  <  4  and  i  =  1,  0  <  n  <  4,  respectively.  It  may  be  noted 
that  these  joint  probabilities  are  all  unimodal  for  t  =  0.  For  :  =  1,  P\I(t)  =  1,  A T(f)  =  0 
approaches  one  monotonically  a s  t  —*  oc.  We  note  that  this  joint  probability  is  also  the 
distribution  function  5iV(o),o(0  °f  t'ime  r(‘(iuire^  f°r  completely  eliminating  all  errors 
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in  the  software  system,  as  shovvm  in  (2. 4).  After  certain  initial  period,  ail  other  joint 
probabilities  with  i  —  1  decrease  monotonicaily  to  zero.  Figure  4.8  depicts  P\I{t)  —  t'j 
for  t  =  0,1.  This  probability  with  i =  1  is  uninodal  having  the  peak  approximately  at 
time  t  —  2  and  decreases  to  zero  as  t  — »  oc.  Correspondingly,  P[I{t)  =  lj  has  minimum 
around  t  =  2  and  goes  to  one  as  t  — ♦  oc. 

We  next  turn  our  attention  to  the  software  system  availability /relia.bility  measures 
based  on  the  past  observation.  We  assume  that  at  time  t  =  20  a  repair  for  the  tenth 
software  failure  is  completed.  In  other  words  we  choose  m  —  10  and  t o  =  20  in  the 
condition  (CD)  given  in  Section  2.  The  associated  probability  vector  J,,  of  (2.5)  can 
be  calculated  using  the  Laguerre  coefficient  vectors  of  /^,.(0-,t)  as 
(4.18) 

ST  =(0.4185245993,  0.24850886S8,  0.1781236373.  0. 104055  v505.  0.050720 

The  counterparts  of  Figure  4.6  and  4.7  are  given  in  Figure  4.9  and  t  10.  •*.. 

joint  probability  of  the  system  availability,  unabailabilitt  ami  the  r  ;n 

software  system  at  time  tu  4-  r  given  (CD).  It  should  be  noted  that  :.*.••  <  ur-.  * 

in  Figure  4.10  is  the  distribution  function  G ,%•(,  «  ; ( r )  of  time  until  co::. ;>;*•:••  tie!.ugg.:,„ 
under  (CD)  given  in  (2.9).  The  joint  probability  of  I{t  )  and  M  {{ )  given  i  CD  >  , i 

in  Figure  4.11,  corresponding  to  the  formula  in  (2.8)  with  ?  -  1  and  ?  5  Fi:i.>.i;%  tin 

survival  function  of  the  time  until  the  next  software  failure  dcn.>tod  bv  11  .  ...  in  (2  10 
is  exhibited  in  Figure  4.12.  We  note  that  lim,--.^  lV(ii  rn(r )  ._  JLI  (  ^  -  0.41852  It  ', 
i.e.  with  this  probability  all  errors  in  the  software  system  have  been  eliminated  by  tim." 
t  -  20  and  there  will  be  no  software  failure. 

Ali  computations  were  done  in  DFC20  in  a  time  sharing  mode  using  APL  as  a  pro- 
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gramming  language.  The  DEC20  APL  implementation  is  the  double  precision  system 
which  uses  a  precision  of  18  decimal  digits.  Relevant  formulas  were  usually  coded  in  a 
straightforward  way  with  no  attempt  made  to  optimize  the  subroutines  for  speed  and 
accuracy.  Evaluation  of  all  performance  measures  presented  in  this  section  required 
approximately  several  minutes  of  CPU  time.  No  evidence  of  numerical  problems  were 
observed  through  the  entire  procedure. 
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Table  4.1  Moment  of  f,  m  n(0~,  t),  i  =  1,  m  =  10 


n 

0th  moment 

1st  moment 

2nd  moment 

0 

0.00356S327439606S2 

0.00356832743960683 

0.0757263591 74633S5 
0.07572635917463717 

1.68068784447817074 
1.68068784447949448  ! 

1 

0.00210648615538208 

0.00210648615538209 

0.04302096760275793 

0.04302096760275967 

0.92238993277688258  ; 

0.92238993277757633 

2 

0.00153425632614903 

0.00153425632614904 

0.0306 19960S3S 17533 
0.03061996083817608 

0.64343981552739176  1 

0.64343981552769234  1 

0.00091263323266573 

0.00091263323266573 

0.01777617904156510 

0.01777617904156599 

0.36542344506696303 

0.36542344506731721 

4 

0.00044773521740952 

J. 00044773521740952 

0.0086S505502745937 

0.0086S505502745974 

0.17791002799323507  j 

0.17791002799338249  i 

upperline:  values  via  recursive  formula 
lowerline:  values  via  Laguerre  sharp  coefficients 


Table  4.2  Moment  of  5,  ^  M  —  41 


n 

:  0th  moment 

1st  moment  >  2  i 

0 

j  5.2021934394 
5.2021934383 

50.9448926126 
50.9448921478  ! 

1 

-1.4812521192 

-1.4812521199 

-15.4420556088 

-15.4420558898 

2 

;  -0.9621601960 
-0.9621601971  i 

;  -8. 13966S3S07  j 

-8.1396688012  i 

3  | 

! 

-0.502233S213 

-0.502233S215 

Lu 

-0.1756699416 
-0.1756699420  ! 

-0.7508504584 

-0.7508505963 

upperline:  values  via  recursive  formula 
lowerline:  values  via  Laguerre  sharp  coefficient 


Table  4.3  1  -  vM {t ),  M  =  4 1 


time 

1  -  t’M (f  ) 

0 

0 

4 

6.505  x  ID-1'1 

8 

6.505  x  10~19 

12 

6.505  x  10“19 

16 

8.674  x  10“19 

20 

6.505  x  10-19 
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NUMBER  OF  FAILURES 
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Appendix 
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In  this  appendix,  we  provide  a  concise  summary  of  the  matrix  Laguerre  transform. 
The  reader  is  refered  to  Sumita(l9S4)  for  more  detailed  discussions.  The  Laguerre  poly¬ 
nomial  Ln{x)  of  degree  n  is  defined  by  the  Rodrigues  formula 


The  corresponding  Laguerre  functions  £n(x)  —  e  -zLn(x),  n  —  0, 1,2, ... ,  for  an  orthog¬ 
onal  basis  of  L2(0,oc)  =  {/  :  —  R  :  /0°°  /: (x)<fx  <  oc}.  The  Lapace  transform  of 


A,(r)  is  given  by 

(A2)  A„(s)  =  f  e~,zCn{x)dx  =  — ^ 

*  0  5  —  0 


s  — 


f  ,  n  =  0,1,2,...,  Ae(s)  > 


We  define  the  linear  space  L„  of  A'  x  A'  matrix  functions  by 


(A3)  L2  -  (g(x)  =  (at;(x))  ;  at;(x)  £  L2  for  all  1  <  t,j  <  A'}. 

It  can  be  readily  seen  that  the  matrix  Laguerre  functions  £n(x)  =  £n(z)  Z  provides  an 
orthonormal  basis  of  I,  where  £  is  a  A'  <  K  identity  matrix.  Then  for  any  g[x)  £  one 
has  the  Fourier-Laguerre  series  expansion 

00  a  4.  roc 

(A4)  g(z)  =  XI  EX!*):  gn  =  j  a{x)[Jx)dx. 

n~‘o  L 

The  pointwise  convergence  of  (A4)  can  be  assured  under  certain  conditions  regarding  the 
smoothness  and  the  rapidly  decreasing  property  of  q(x).  Let 

4.  4.  4. 

(A5)  =  5o  ’  £«=E„-£n-l’ 

4. 

and  define  the  two  matrix  generating  functions  A  =  y^L0a^un  and  7*  -  E“=0s£u\ 
We  note  that 

(A6)  ?T(u)  =  (i  -  ti)7i(u),  |uj  <  1, 

iii 


7T 


—  „  -  ,  -  «  -  v  ■* y  *  jf  p  v  ■  v  ■  y  w  y  m [V  •  l  *  V  v  9  v.1  *  V'  w  V  v  IT*  7* »)■  ii  m  v  rnr rt  7 ■v  »_nf  ^ t 


f* 


pa 


K-‘- 


fcv 


r*  -■ 


tVv 

I 


i 

*  -■■  .' 

k\-.V 

K--'V 

>  .■-  ,■ 

[••-V 

■  •. 


W-.  ■-- 

:■/; 


and 

(-47)  in=  fl  %■  ri=0'1’2 . 

m  =  u 

From  (A2)  and  (A-4),  one  sees  that 

(.48)  a(s)  =  [  e~!Za(x)dz  —  YJ  q_n  - 

Jo  -s 


a.  1 


By  letting  u  = 


(AO) 


n=0 

this  leads  to 


i  U  -  k 


Re(s)  >  — . 


7r(«)  =  s(J 1  “ 


>2  1  -  u. 

Equation  (A9)  is  the  key  formula  for  the  matrix  Laguerre  transform,  providing  a  bridge 
between  continuum  operations  and  lattice  operations.  For  the  matrix  convolution  c(x)  = 
Jo  a(x  -  y)b{y)dy  with  q(x),  £(x)  £  for  example,  it  can  be  readily  seen  that 


MlO) 

or  equivalently, 
(.411) 


tc-(u)  =  tz(u)t;m. 


c*  =  4*. 

=n  —  J  —] 

;=o 


The  matrix  Laguerre  transform  maps  matrix  functions  a(x),  £(x)  £  L ;  into  the  matrix 
sequence  (a“)o°  and  (4^)o°-  Correspondingly  the  matrix  convolution  on  continuum  is 
mapped  into  the  lattice  convolution.  The  resulting  sharp  matrices  (g* )°°  can  be  converted 
to  the  Laguerre  dagger  coefficient  matrices  (m  using  (AT).  The  values  £n(x)  can  be 
generated  efficiently  via  the  recursive  formula 

(A12)  fn+j(x)  =  — ~:(2n  1  -  x)f„(x)  -  n£tl_l{x).,  n>l 

v  n  —  1 

j  1 

starting  with  £q(x)  =  e~  i1  and  fi  (x)  =  (1  -  x)c“:1.  Hence  the  Laguerre  dagger  coefficient 
matrices  (c  1)3°  can  be  inverted  back  onto  continuum  via  the  series  representation  (A4). 
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Other  continuum  operations  arc  mapped  into  lattice  operation  in  a  similar  manner.  We 
list  only  a  few  operational  properties  below. 

Integration 

Let  £(x)  =  a(y)dy.  Then 


{A  13) 


-Z  a 
=n 


Multiplication  bv  polvnomir.als 
Let  g(z)  =  xg(x).  Tiien 

(AM)  <T  =  -A:;(n  -  n  >  0. 

where  A  a  —a  —a  , . 

=n  =rn-l 

The  first  two  moments  of  a(x)  can  be  also  evaluated  in  terms  of  (a“)£L0. 

Moment  formula 

-DC  OC 

(A15)  /  zla(z)dz  =  4*  V]  ( - ,  0  <  t  <  2. 

0  n=0 

The  higher  moments  can  be  comuputed  by  combining  ( A 1 4 )  and  (A15)  repeatedly.  The 
moment  formula  (A15)  provides  a  empericai  tool  for  deciding  the  truncation  point  of 
the  series  representation  (A4),  when  the  moment  values  are  known.  Extensive  numerical 
experiments  suggest  that  if  the  truncation  point  is  chosen  to  satisfy  a  given  accuracy  for 
the  moment  value  in  (A  15),  then  the  series  representation  ( A 4 )  also  satisfies  the  same 
accuracy. 

It  should  be  noted  that  for  many  of  matrix  functions  of  interest  in  applied  probabihty 
and  statistics,  the  Laguerre  coefficient  matrices  can  be  obtained  either  numerically  or  via 
certain  numerical  procedure  based  on  (A9). 
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